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Section 0. Introduction 


Methods for modeling high speed propulsion systems will be discussed. Included in 
this category are internal flow propulsion systems that do not contain rotating machinery. 
Specifically, inlets, ramjets, and scramjets are considered. Some direct extensions to rotating 
machinery are available, however they are not discussed. It should be stressed here, that the 
main application of the resulting models is some aspect of control system design. This then 
requires a clear understanding of the trade-offs between model complexity, to correctly model 
the system, and model simplicity, to allow control system design. This report then clarifies the 
modeling assumptions that are necessary for various levels of model complexity. 

Gasdynamic flows are usually described by the set of Euler equations (see the next 
page or [Hirsch; Anderson, Tannehill, & Pletcher]). As internal flows are the primary 
consideration, a quasi-one-dimensional approximation is often satisfactory and is pursued here. 
Very high speed flows sometimes require the inclusion of viscosity effects, temperature 
dependent gases, and chemical reactions. Temperature dependent gases are easily accounted 
for in the quasi-one-dimensional Euler equations by allowing the particular ratio of specific heat 
to be a function of temperature [Anderson; Ames Tables]. Viscosity effects require the addition 
of an extra term to the momentum equation of the Euler equations to change them into the 
Navier-Stokes equations [Anderson, Tannehill, & Pletcher]. This extra term is easily 
accounted for by all of the methods discussed, as long as two dimensional effects such as 
boundary layers and turbulence are not considered [Lin]. Chemically reacting gases are more 
of a problem, as each chemical species requires the addition of another continuity equation, as 
well as diffusion, energy release, and reaction rate terms to the Euler equations [Oran & Boris]. 
Furthermore, [Anderson] shows that hydrogen in air requires seven species and possibly 140 
reactions. Although this is important, the increased complexity will not be considered practical 
at this time. Chemical reactions can be accounted for by the addition of heat in the energy 
equation for the present. Clearly the addition of all these effects into either a nonlinear or 
linearized model is an important area for future research. Other important phenomena that will 
not be discussed here are boundary layers, turbulence, oblique shocks, inlet buzz, and 
combustion instabilities. These phenomena require at least two space dimensions for an 
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accurate and physically realistic representation [Hsieh, Wardlaw, & Collins; Lin]. Extending 
the discussion contained here to these higher dimensional phenomena is of great future 
importance. 

The quasi-one-dimensional Euler equations are [Varner et al (LAPIN Report)] 


continuity 


d(pA) , 5 (P uA ) _ M 

8t dx 


momentum 


d(puA) + 3[A(p+pu 2 )] _ p 3A + p 
3t 3x dx 


3(EA) 3[Au(E+p)] 3A . 

energy — — - + — — - — — = -p — + Q 

3t dx dt 

where p = pRT = (y-l)[E - 0.5pu 2 ]; E = p[c v T + 0.5u 2 ]; y = ~ ; and the other variables 
have their usual definitions [Varner et al (LAPIN Report)]. The solution of these equations 
with various levels of complexity is the purpose of the following discussion. As mentioned 
earlier, many high speed propulsion systems can be approximately represented by these 
equations. 


The discussion which follows is separated into four areas. These are in order of 
discussion; 1) CFD models that represent the entire nonlinear system, or high order nonlinear 
models; 2) high order linearized models derived from the fundamental physics; 3) low order 
linear models obtained from the other high order models; and 4) low order nonlinear models 
(order here refers to the number of dynamic states). Included in the discussions on modeling 
will be any special considerations based on the relevant control system designs. Where 
necessary, some digression into specialized control techniques will also be undertaken. 
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This section considers nonlinear models that contain large numbers of states and are 
often very accurate. Most of these methods are based on computational fluid dynamics (CFD) 
and are typically finite difference representations of the quasi-one-dimensional Euler equations 
above. These methods are further divided into methods of high, medium, and low accuracy. 

High Accuracy Methods 

The methods considered to be highly accurate are usually implicit methods [Varner et al 
(LAPIN Report); Hirsch; Anderson et al]. They consider more than one space dimension and 
can have millions of spatial lumps, and thus millions of states. Although these methods are 
perhaps not practical for control system design, or real-time simulation, they do provide a very 
good representation of the system dynamics. Currently, methods of this type will require 
hours, or even days, of computing time to provide a complete transient response. In the near 
future, however, increased computational speed could allow real-time CFD models to be 
running in parallel with actual high speed propulsion systems. With appropriately placed 
sensors, the real-time CFD model would function as an observer, providing a control algorithm 
with information about all the flow parameters everywhere in space. This will then push the 
control system designer to find a way of incorporating all this information into a control 
system. The question being, "if we have all this information about all the states of the flow, 
everywhere, what then do we do with it in order to provide an improved control system 
design." This question has been considered in [Hartley & Xia] where distributed functional 
observers and controllers are shown to obey a separation principle. This problem will also be 
discussed in Section 2. 

! Medium Accuracy Methods 

! 

I Medium accuracy methods are in many ways more practical than the high accuracy 

j 

methods discussed above. They use the quasi-one-dimensional Euler equations and finite 

j difference schemes that have a local truncation error proportional to the timestep squared and/or 

the spacestep squared, and are thus second order accurate. These methods typically contain 

| 
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hundreds of states and are thus somewhat more tractable and manageable. The method of 
choice in this category is clearly MacCormack's method [Hirsch; Anderson, Tannehill, & 
Pletcher; Hartley, Melcher, & Bruton], or any of the other variations on the two-step Lax- 
Wendroff idea. This method uses a predictor-corrector approach where an approximate 
solution is predicted in a downstream fashion, and the improved solution is corrected in an 
upstream fashion. The reader is reminded that the main difficulty in simulating the quasi-one- 
dimensional Euler equations is the first spatial derivative appearing there. The finite difference 
replacement must allow information to travel in both directions. The naive choice of either 
forward or backward differences does not allow this simultaneously. The next obvious idea of 
using central differences is unfortunately unstable. Thus one is left to be clever. The medium 
accuracy clever solution is the family of two step Lax-Wendroff methods [Hirsch; Peyret & 
Taylor]. 

MacCormack's method has been very useful in modeling the Euler equations. It has 
been used to model various propulsion system components [Varner et al (LAPIN Report); 
Hartley, Melcher, & Bruton]. The user is reminded that smoothing is required after the 
predict-correct mode in order to reduce the spatial oscillations due to the shock discontinuity 
[Interim report]. 

Low Accuracy Methods 

Here low accuracy will refer to methods that have a local truncation error on the order 
of the timestep or spacestep size, or first order accurate methods. The most common methods 
in this category are the flux splitting methods [Anderson, Tannehill, & Pletcher; Hirsch]. 
There are several variations of these methods, based on the particular local splitting method, 
but the general idea is that somehow the flow at a given lump is separated into that which is 
traveling upstream and that which is traveling downstream. Once this is done, backward 
differences can be used on that which is traveling downstream and forward differences can be 
used on that which is traveling upstream. This then preserves accuracy, stability, and shock 
capturing while still only using first differences. In real-time simulation, these splitting 
methods are no faster than MacCormack's method, due to the necessity to separate the flows at 
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each lump on each timestep. Even with pre-separated methods, multiple function evaluations 
are still necessary and thus little speed is gained. 

An alternative to these splitting methods is termed physical lumping. Here, the spatial 
differencing of each term in the given derivative is done in a physically intuitive sense. This is 
discussed in more detail in [Immel, Hartley, & DeAbreu]. Although mentioned in [Roache] 
and attributed to [Courant, Isaacson, & Rees], the approach has received little attention in the 
CFD community. This is probably due to the difficulty associated with predicting accuracy and 
stability. However [Immel, Hartley, & DeAbreu] have shown that reasonable accuracy and 
shock capturing are possible in a supersonic inlet model. More study should be given to this 
approach as it provides a relatively simple and inexpensive method for simulating the quasi- 
one-dimensional Euler equations. Note that both the splitting methods and the physical 
lumping methods also require hundreds of states, as did the medium accuracy methods. 


Section 2. High Order Linear Models 

Where the accurate large perturbation models of the last section are important for 
control evaluation; small perturbation, or linearized, models are usually required for control 
design. The next two sections of this document will address techniques for obtaining 
linearized models of the given Euler equation flowfield. This section will address techniques 
for obtaining linear models with a large number of dynamic states, usually based on physical 
first principles. Section 3 will address how to reduce the size of these to something that is 
more reasonable for control design. The three basic methods discussed in this section are the 
approximate linear solution of the exact nonlinear system, exact solution of an approximate 
linear system, and approximate solution of an approximate linear system. All of these methods 
first require the particular system to approach the desired steady state. Once the steady state 
solution is obtained, a linear model which is only valid for small perturbations from this 
operating point can be created. 
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CFD Based Methods 


CFD based methods provide good nonlinear models for the Euler equations. 
Fortunately, the flux splitting methods also allow a simple linearization technique to be used to 
create a high order linear approximation to the Euler equations. This method again requires the 
steady state conditions of the flowfield. Once these are obtained, either computationally or 
analytically from the design specifications, the desired number of spatial lumps must be 
determined. Usually the number of spatial lumps is determined from the desired spatial 
accuracy and the length of the flowfield. Given this, a particular CFD method is selected. The 
flux splitting methods are particularly suited to this approach. Mathematically combining the 
CFD method with the Euler equations and linearizing, allows the construction of a block 
tridiagonal system matrix, with the number of states equal to three times the number of lumps. 
Clearly, to obtain reasonable spatial accuracy requires a large number of lumps, perhaps 
hundreds. This approach, using one of the split flux methods, is discussed more completely in 
[Chicatelli; Chicatelli & Hartley]. It is discussed using the physical lumping method in [Immel, 
Hartley, & DeAbreu]. 

Exact Solution of the Linear Flowfield 

Since most propulsion systems have variable cross-sections, the properties of the flow 
depend upon spatial position. A linear approximation to the Euler equations themselves then 
requires some additional assumptions. Either many linear lumped volumes must be used and 
pasted together, or a single linear lumped volume can be used with volume averaged flow 
properties, or the linear flowfield must be solved with spatially dependent flow properties. The 
latter is extremely difficult to do generally and is not considered further. The second must be 
done before the first, so the second will be considered first in this subsection. This is 
essentially obtaining an exact solution to an approximate problem. 

To use this approach, the volume averaged steady state solution must be found for the 
entire flowfield, which is considered to be one lump. Then, the linearized flowfield is 
described by a linearized set of the Euler equations, with the appropriate boundary conditions. 
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inputs, outputs, and shock discontinuities also linearized. The resulting system can then be 
Laplace transformed and the resulting spatial boundary value problem can be solved, letting the 
Laplace s-variable be a constant. The end product of this procedure is then a non-rational 
meromorphic transfer function (or transfer matrix) which basically has an infinite number of 
states. When done correctly, this transfer function represents the Laplace transform of the 
Green's function of the flowfield and is very useful for studying the input-output properties of 
the system [see Butkovskii for more general information on this approach]. This approach has 
been effectively applied to the isentropic subsonic flow problem (a modeling assumption that 
eliminates the energy, or entropy, equation in the Euler equations) with the results presented in 
[Sarantopoulos & Hartley]. This same procedure is currently being expanded to consider 
supersonic- subsonic combination flows, supersonic flows with detonation, and multiple lumps 
pasted together [Sarantopoulos]. 

Approximate Solution to the Linear Flowfield 

This approach basically follows that above, but the meromorphic transfer functions are 
not obtained in the process. Before the algebra is performed, the flow is diagonalized, as in 
flux splitting, and the resulting conservative information delays (transport lags) are replaced by 
Pade approximations. This is the method of [Willoh; Cole & Willoh] and has been used 
extensively. Since the meromorphic time delays have been replaced by finite order Pade 
approximations, the resulting number of states is equal to 3*(number of lumps)*(Pade 
order)+l. Another method that has been used, and is apparently more accurate and 
appropriate, is to replace the time delays by discrete delays (z-inverses) in a digital simulation 
of the entire system [Hartley]. Note that the general Cole-Willoh approach does not contain the 
spatial input-output information of the Sarantopoulos approach above. Although certainly 
useful, it would probably be more accurate to approximate the entire Sarantopoulos transfer 
function than to approximate little pieces of the entire transfer function and then multiply them 
together as in the Cole-Willoh approach. 
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Based on the methods of the previous section, several model reduction approaches are 
possible. These are discussed first in this section. Following this, some other approaches for 
generating low order linear models are presented. 

Reduction of Linearized CFD Models 

The linearized state space models based on the CFD approach discussed in the last section are 
very well suited to the wide variety of modern linear model reduction techniques based on 
balancing [Moore; Laub, Heath, Paige, & Ward; Safanov & Chiang], These techniques have 
been applied to both the split flux and the physical lumping approaches with considerable 
success using the MATLAB 'schmr' function. An 80-90% reduction of states has been 
possible in supersonic inlet models as well as ramjets models with upstream, downstream, and 
mid-stream inputs [Chicatelli; Chicatelli & Hartley; Immel, Hartley, & DeAbreu]. Essentially, 
many states are only slightly observable/controllable. This is particularly true for the 
supersonic part of the flow with an input in the subsonic region. Clearly, the states in the 
supersonic region should not be controllable from downstream. Furthermore, there is some 
tendency to lose controllability/observability based on the particular positions chosen for the 
actuators and the sensors. These isolated state problems must be more clearly understood 
before real-time CFD observers can be effectively implemented. To this end, an input-output 
separation principle for observing and controlling spatially distributed systems has been 
developed and studied in [Hartley & Xia]. 

Direct Reduction of Meromorphic Transfer Functions 

The idea here is to obtain some finite order system based on the Sarantopoulos-Green's 
transfer function approach discussed in the last major section. This was already introduced via 
the Cole-Willoh approach where the lump time delays were replaced by either Pade 
approximations or discrete-time delays. The idea here will be to replace the entire transfer 
function by an approximation. Several methods are available for approximating general 
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meromorphic transfer functions. Some of the more common methods are truncated Taylor 
series approximations of the entire transfer function (effectively a Pade approximation of the 
meromorphic function and not necessarily stable [Edrei, Saff, & Varga]), truncated Taylor 
series expansions of the numerator and denominator separately (only stable if done carefully 
with positive exponential powers rather than negative powers), infinite product expansions (if 
the transfer function is a simple function), and partial fraction expansions (only if enough 
approximate pole positions can be found). More information on all these methods can be 
found in [Henrici]. An alternative to this approach that has only recently been developed is 
based on approximation of the system Hankel operator [Partington]. This method is generally 
very difficult and requires knowledge of the system impulse response. An approximation of 
this method has been applied to representative inlet systems in [Hartley & DeAbreu] using a 
numerical inverse Laplace transform. Another approximation has been presented by [Gu & 
Khargonekhar, & Lee] which directly uses an inverse Fourier transform. 

Optimized Low Order Models 

The idea here is to optimally choose the parameters of a given reduced order model by 
reducing some performance measure. This is fairly simple conceptually, and many useful 
performance measures are available. If a time response is available, either from the transfer 
function or from time response data, a time domain performance measure, such as squared 
error or time weighted squared error can be used. This can be also considered an input-output 
identification method. An exhaustive study using this approach can be found in [Piercy] where 
it is shown that the batch total least squares method is somewhat more robust and accurate than 
the usual batch least squares method when applied to a supersonic inlet problem. Alternatively, 
if only the transfer function is available, minimization of some performance measure in the 
frequency domain is possible. The reader is reminded that this is not usually a linear problem 
and can have many local minima. This approach has been used to reduce the order of a 
representative inlet system in [Dariush & Hartley]. 

Other Methods for Generating Low Order Models 

Obviously other methods are available for generating linear low order models of high 
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speed propulsion systems. These are reviewed more completely in [Hartley], One approach 
worth mentioning is that of representing the Euler equations by several lumped circuit 
equivalents [Stalzer & Fiedler]. If an isentropic equivalent is used, the usual transmission line 
equations result. Otherwise another parallel voltage path is required for the third Euler 
equation. This circuit approach allows for quick analysis and testing on the currently available 
circuit analysis packages, such as SPICE. 


Section 4. Low Order Nonlinear Models 

Although there are many model reduction methods available for linear systems, 
unfortunately there are very few available for nonlinear systems. Hence it is generally difficult 
to start with the CFD models of Section 1 and obtain reasonably accurate, low order, large 
perturbation nonlinear models. This method is somewhat simplified for low speed subsonic 
flow where a variety of physical lumping methods have been used [Hartley; Krosel & Bruton; 
Colboume]. It is also fairly simple for high speed supersonic flow [Interim Report: Chicatelli 
Scramjet]. Unfortunately, the inclusion of shock dynamics makes this problem much more 
difficult. The problem with the usual model reduction methods, and low order physical 
lumping methods, is that the states in a nonlinear finite difference CFD model are not always 
readily available in a closed form (such as MacCormack's method), and when they are (as in 
physical lumping) they each represent information in a given point in space. When the shock 
undergoes a large perturbation, it is not intuitively clear to individuals, or the reduction method, 
which combinations of states are important. Currently, research is continuing in this effort 
using the physical lumping models as they allow easy access to the states. 

The creation of single closed form low order nonlinear models which include the entire 
range of dynamics of an inlet, ramjet, or scramjet, has not yet been accomplished. First, the 
model must be capable of transitioning from one type of behavior to another, and allowing both 
behaviors to simultaneously coexist. For example, the shock in a supersonic inlet can be in a 
started mode, which corresponds to a point attractor for the state space system, or it can be in 
buzz mode, which corresponds to a limit cycle for the state space system. Usually both 
behaviors are possible simultaneously; the mode being dependent upon the initial conditions 
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and the magnitude of any perturbations. The preservation of multiple operating 
conditions/attractors and the maintenance of their dimension is addressed further in 
[Mossayebi, Hartley, & DeAbreu-Garcia; and Hartley, Killory, DeAbreu-Garcia, Abu 
Khamseh]. 

Four techniques are suggested for future consideration. One, the method of [Martin] 
allows transitioning through preprogrammed logic. When a transition is called for, a new 
model is inserted for the old one. Although crude, the method certainly works. Second, 
numerical optimization of a low order nonlinear model from very accurate CFD models of all 
necessary phenomena is considered. Here, it would be necessary to have a two or three space 
dimensional model to generate the data, as most instabilities fundamentally require two or more 
space dimensions to represent them. Also, it would be necessary to chose an appropriate 
model structure which would not be readily apparent initially. This is probably an approach 
where neural networks could work quite well due to the large amount of uncertainty. Third, a 
model reduction technique that appears to be applicable to nonlinear systems is aggregation 
[Aoki], This method essentially allows the user to choose the combination of any desired 
states to keep for the reduced order model. It is then possible to algebraically eliminate the 
remaining combination of states. This research is currently being pursued. Fourth, spectral 
methods for the Euler equations have a tremendous potential for generating these accurate low 
order large perturbation models [McCaughan; Culick, Lin, Jahnke, & Sterling]. These 
methods assume an infinite orthogonal set of spatial eigenfunctions for the spatial boundary 
value problem. Using separation of variables and some cancellation, the temporal coefficients 
for the spatial eigenfunctions form a low order set of nonlinear ordinary differential equations 
for a subset of the spatial eigenfunctions. The resulting simulation then would yield a spatially 
continuous solution which would change as the temporal coefficients changed. The major 
drawback of this approach is that, as usual, hyperbolic systems like the Euler equations do not 
readily lend themselves to this approach. More information can be found in [Gottlieb & 
Orszag]. It should be noted that this approach has been very effective for other flow systems 
including turbulence [Berge, Pomeau, Vidal; McCaughan]. 
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Section 5. Conclusions 


A unified systematic modeling approach has been presented for the modeling of high 
speed propulsion systems. When choosing a particular modeling approach, the decision 
between phenominalogical accuracy and applicable complexity must always be made. This is 
true for both control design and for control evaluation, or real-time simulation. The methods 
discussed here are for the quasi-one-dimensional Euler equations of gasdynamic flow. The 
basic methodology and organization applies, however, to any other nonlinear spatially 
distributed system. The essential nonlinear features accurately represented by the quasi-one- 
dimensional Euler equations and the modeling methods discussed here are large amplitude 
nonlinear waves, including moving normal shocks, hammershocks, simple subsonic 
combustion via heat addition, temperature dependent gases, detonations, and thermal choking. 
For accurate representations of oblique shocks, inlet buzz, boundary layers, flow separation, 
turbulence, and combustion instabilities, it would be necessaiy to consider at least two space 
dimensions in the Euler equations and viscosity would need to be added. This is the next 
logical step in the development of modeling methods for the control of high speed propulsion 
systems. 
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